Workshop 

Future  Directions  in  Systems  and 
Control  Theory 

Thursday ,  June  24 
Viewgraphs  -  Volume  3 


DISTRIBUTION  STATEMENT  A 

Approved  for  Public  Release  | 
Distribution  Unlimited 


20000118 125 


jfEOQP&MFi  nfSeaCEBD  1 


A  Framework  for  Control, 
State  Estimation,  and 
Verification  of  Hybrid  Systems 


Manfred  Morari 

Alberto  Bemporad 
Giancarlo  Ferrari-Trecate 
Domenico  Mignone 
Fabio  Danilo  Torrisi 


Automatic  Control  Laboratory 
Swiss  Federal  Institute  of  Technology  (ETH) 

moraxiSaut . ee . ethz . ch 


Technical  Reports:  http://ww.aut.ee.ethz.ch 


mmmmm  m  i 

gm  g  gg  Miboti  NHAMArtj 


Premise 


•  The  discussed  problems  are  important  but 
inherently  difficult 

•  All  useful  techniques  must  involve  significant 
off-line  and/or  on-line  computation 

Results 

•  New  System  Type: 

Mixed  Logical  Dynamical  (MLD)  System 

•  Many  practical  problems  can  be  represented 
in  MLD  form 

•  Control,  estimation,  and  verification 
require  solution  of  Mixed-Integer  Linear  (or 
Quadratic)  Programs  (MILP,  MIQP)  for 
which  efficient  techniques  are  becoming 
available 
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Example:  Verification  of  an 
Automotive  Active  Leveler 
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•  For  a  given  set  of  initial  conditions,  variation  of 
parameters,  and  road  conditions: 

(VP#1)  Find  worst  performance 

(VP#2)  Verify  that  height  is  within  desired  limits 
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Example:  Control  of  a 
Hydroelectric  Power  Plant 


Objective: 

•  Maximize  power  generation 

Manipulated  variables: 

•  Turbine  ON  /  OFF 

•  Stepper  motors  (-1,0,1) 
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Hydroelectric  Power  Plant 


Constraints 

•  on  actuators  (flaps,  gates) 

•  on  level 

Targets  (  =  end  point  constraints) 

•  level 

•  actuators  (flaps,  gates) 

Preference  hierarchy  /  heuristics 

1.  Flaps  rather  than  gates 

2.  Turbine  ON  >  Tm*n 

3.  Stepper  motor  ON  >  Tmin 

4.  Transitions:  turbines  /  flaps  /  gates 

5.  Stepper  motor  OFF  <  Tmax 
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Example:  Three  Tank  System 


Nominal  behaviour: 

•  Liquid  level  in  tank  1  controlled  by  pump  1 

•  Liquid  level  in  tank  3  controlled  by  switching 
valve  V\ 

Faults: 

Leak  in  tank  1 
4>2‘  Valve  Vi  blocked  closed 

Maintaining  nominal  operating  condition 
requires  fault  detection  &i  control 
reconfiguration 
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Example:  Gas  Supply  System 

(Kawasaki  Steel) 


•  Maintain  hold-up 

•  Maximize  use  of  gas  for  combustion 

•  Limitations  on  gas  flows,  heat  produced,  gas 
volume  in  holders 

•  Complex  rules  for  boiler  switching 
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Premise 


•  The  discussed  problems  are  important  but 
inherently  difficult 

•  All  useful  techniques  must  involve  significant 
off-line  and/or  on-line  computation 


Results 

•  New  System  Type: 

Mixed  Logical  Dynamical  (MLD)  System 

•  Many  practical  problems  can  be  represented 
in  MLD  form 

•  Control,  estimation,  and  verification 
require  solution  of  Mixed-Integer  Linear  (or 
Quadratic)  Programs  (MILP,  MIQP)  for 
which  efficient  techniques  are  becoming 
available 
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General  Class:  MLD  Systems 


General  formulation  of  Mixed  Logical-Dynamic 
(MLD)  systems 


x(t  +1)  =  Ax(t )  +  Biu(t )  +  B28{t)  +  -632(f) 
y(t )  =  Cx(t )  +  Di«(f)  +  D2S(t)  +  632(f) 
E25(t )  +  632(f)  <  £7x^(t)  +  64a  (t)  +  £5 


•  State  : 


*  =  '  ice  Rnc,  xt  e  {0, l}n* 

Output  : 

V  =  M  ’  Vc  G  RPc,  y^  e  {0, 1  yt 

Vi_ 

Input  : 

u=  Uc  ,  uc€  Rmc,  ut  E  {0,  l}me 
uej 

Auxiliary  binary  variables  :  <5  €  {0,  1)6 
Auxiliary  continuous  variables  :  2  €  Mrc 
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Assumptions 


Recall:  x(t  +  1)  =  Ax(t)  +  .Biu(£)  +  jE^W  -\-Bsz(t) 

y(t)  =  Cx(t)  +  Diu(t)  +  D2S(t)  +  D3z(t) 
E25(t)  +  E^zit)  <  Eiu(t)  +  E4x(t)  +  E5 

•  Wei  I- Posed  ness  : 

{a:(t),  u(t)}  — >■  {x(t  +  l)}  single  valued 
{x(t),  u(t)}  -»■  {y(t)}  single  valued 

This  allows  to  define  trajectories  in  z-space 

•  Physical  Constraints  : 

€  Mnc+mc  :  Fxc  +  Guc  <  j 

C  is  a  bounded  polyhedral  set  (not  restrictive 
in  practice) 

This  allows  to  define  upper  and  lower  bounds 
M ,  m 

Note:  Physical  constraints  are  included  in 
the  inequality 
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Propositional  Calculus  and 
Linear  Integer  Programming 


Truth  table 
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Propositional  logic  Integer  linear  inequalities 
X  e  {F,  T}  S  e  {  0,  l  } 


Xi  v  X2  is  equivalent  to  <5i  +  62  >  1 

X\  a  X2  is  equivalent  to  <5i  +  82  >  2 

~  X\  is  equivalent  to  <  0 
Xi  X2  is  equivalent  to  <5i  —  82  <  0 

X\  X2  is  equivalent  to  8\  —  62  =  0 
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Logic  Facts  Involving 
Continuous  Variables 


[/(a:)  <  0]  A  X 

I 

rH 

+ 

\  1 

1 

VI 

1 

r-s 

B 

[f(x)  <  0]  V  X 

f(x)  <  MS 

~  [/(^)  <  0] 

f(x)  >  € 

[/( *)<0]  ^  x 

— /(as)  +  e  <  (e  -  m)8 

[/(a:)  <  0]  ++  X 

J  f(x)  <  M(l  — 5) 

1  /(a:)  >  e  +  (m-e)5 

•  M,  m  =  upper-,  lower-bounds  on  f(x) 

•  e  =  small  tolerance  (e.g.  machine  precision) 

•  X  E  {F,  T},  xeW1,  tie  {0, 1} 

if /(x)  linear  =>  Mixed-Integer  Linear  Inequalities 
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Product  of  variables 


S3  =  5\52 

-5i  +  63  <  0 

—82  83  <  0 

81  +  82  ~  83  <  1 

y  =  5f(x ) 

1 

f  y  <  MS 
\  y  >  mS 

|  y  <  f(pc)  —  m(  1  —  5) 

U  >  /(»)  —  M(1  —  <5) 

•  M,  m  =  upper-,  lower-bounds  on  /(*) 

•  e  =  small  tolerance  (e.g.  machine  precision) 

if /(*)  linear  =*>  Mixed-Integer  Linear  Inequalities 
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An  Example 


System: 

x(t  + 1) 

y(t ) 


0.8 


a(t)  = 


Rewrite  as 


cos  a(t)  —sin  a(t) 
sina(t)  cosa(t) 

[1  0]x(t) 

|  if  [1  0]*(i)  >  0 
|  if  [1  0]x(f)  <  0 


x(t)  + 


u(t) 


x(t  +  !)  =  [/  i]  z(t ) 
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•  The  origin  is  globally  asymptotically  stable 
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System  Theory 
for  MLD  Systems 


•  We  1 1- Posed  ness 

•  Stability 

i: 

!  •  Reachability/Controllability 

f 

j  •  Observability/Reconstructibility 

I 

■  •  Optimal  Control 

I  •  State  Estimation 

I 

J 
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Observability  Test  for 
Autonomous  MLD  Systems 

•  DEFINITION  of  incremental  observability: 

Vxi,  X2  €  T(0) 

T-l 

1 |y(t> *1)  -  yd,  x2)lloo  >  w\\xi  -  T2II1 

t= 0 

•  GOAL  :  compute  minimum  number 
of  output  measurements  necessary  to 
distinguish  two  initial  states  xi,x2  €  Af(0) 

•  TEST  : 

0.  Fix  wm in  as  small  as  desired 

O  <  ■u;rnin=practical  indistinguishability) 

1.  Fix  Tm3X=maximum  number  of  measurements 

(T  >  Tmax=practical  indistinguishability) 

2.  T  <—  1 

3.  Minimize  Yh=o a^i)— ^2) ||oo— ^min ^2|| x 

4.  If  minimum  <  0,  increase  T. 

5.  If  T  >  Tmaxi  STOP  (system  unobservable) 

6.  Go  to  3 

=^>  Mixed-Integer  Linear  Program 
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Observability  Index:  Example 


xi 

x2 
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if  e  <  x\ (t)  <  1,  e  >  0 
otherwise 


y(t)  —  xi  (t), 


e  <  xi(0 )  <  1 


logi 
log  1.1 
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•  Observability  index  T  does  not  depend  on 
state  dimension  n 

•  T  can  be  arbitrarily  large  (also  T  =  oo) 


J 
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System  Theory 

for  MLD  Systems 


•  Well-Posedness 

•  Stability 

•  Reachability/Controllability 

•  Observability/Reconstructibility 

•  Optimal  Control 

•  State  Estimation 
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Model  Predictive  Control 


Reference 


•  MODEL  :  a  model  of  the  plant  is  needed  to 
predict  the  future  behaviour  of  the  plant 

•  PREDICTIVE  :  optimization  is  based  on  the 
predicted  future  evolution  of  the  plant 

•  CONTROL  :  control  complex  constrained 
multivariable  plants 
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Optimization  Problem 


Compute  the  optimal  sequence  of  manipulated 
inputs  which  minimizes 

•  tracking  error  =  output  -  reference 

•  subject  to  constraints  on  inputs  and  outputs 

On-line  optimization  =>-  Receding  Horizon 
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Receding  Horizon 


nast  ▲  future 


•  Optimize  at  time  t  (new  measurements) 

•  Only  apply  the  first  optimal  move  u(t ) 

•  Repeat  the  whole  optimization  at  time  *  + 1 

Advantage  of  on-line  optimization: 

FEEDBACK! 
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Investment  plans 


Predictive  Controller 

•  On-line  optimization  problem  (MIQP) 

T-l 

min  J(vj_1,x(i))  =  T^|K*0  -  Hellg,  +  POI *)  ~  ^||q2 

^  k=  0 

||z(fc|t)  -  Ze\\%3  +  H*(fc|t)  -  Xe\\ Qt  +  ||y(fc|t) 


Ve  lift 


subj.  to  < 


x(T\t)  =  Xe 

x(k  +  l|t)  =  Ax(k\t)  +  Biv(k)+ 

B25{k\t )  +  B3z(k\t) 

y(k\t)  =  Cx(fe|t)  +  Dtv(k)  +  D25(k\t)  +  D3z{k\t) 
E28(k\t)  +  E3z(k\t)  <  Eiv(k)  +  EAx(k\t)  +  E5 


According  to  a  receding  horizon  philosophy, 

u(t)  =  «t(0) 

Repeat  everything  at  time  t  + 1 


Feasibility  of  MIQP  <§t  =  0  =>•  closed-loop 
stability 

Global  optimum  not  needed  ! 
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Closed-Loop  Example 


jc-space 


y(t),  r{t) 


sina(t)  cos  a(t) 

x(t) 

l  if  [1  0]a;(t)  >  0 
\  if  [1  0]x(t)  <  0 
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Computational  Aspects 


•  The  on-line  optimization  problem  is  a 
Mixed-Integer  Quadratic  Program  (MIQP) 

•  No  need  to  reach  global  optimum  (see  proof 
of  the  theorem) 

•  Available  methods:  Generalized  Benders’ 
Decomposition,  Outer  Approximation, 
LP/QP  based  branch  and  bound,  and 
Branch  and  Bound. 

•  Branch  &i  Bound  is  the  most  effective 

•  General  purpose  B&B  MIQP  solvers 
available,  for  both  dense  and  sparse 
problems  (e.g.  Fletcher-Leyffer’s) 
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System  Theory 
for  MLD  Systems 


Well- Posed  ness 
Stability 

Reachability/Controllability 
Observability/Reconstructibility 
Optimal  Control 
State  Estimation 
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Fault  Detection  and  State  Estimation 

for  Hybrid  Systems 

•  Use  Moving  Horizon  Estimation  for  MLD 
Systems  (dual  of  MPC) 

•  Mixed  logic  dynamic  fault  (MLDF)  form 
(Bemporad,  Mignone,  Morari,  1998): 

x(t  1)  =Ax(t)  B\u(t)  B28(0)  “I-  B3z{t) 

+ Bern + m 

y(t )  =Cx(t)  +  Dxu(t)  +  D26(t)  +  D3z{t)  + 

+  D6(/>(t) + at) 

E28(t)  +  E3z(t)  <  ElU(i)  +  E4x(t)  +  E5  +  E6<f>(t) 


Faults:  4>{t)  e  {0, 1}^ 

Disturbances:  £(i)  e  Rn,  C(t)  c  -Rp 


•  At  each  time  t  the  estimates  x(t )  are 

obtained  by  minimizing  an  MIQP  over  a 
horizon  extending  backwards  in  time 


HI/OIO  -3/(fe)llQ„  +  -r'(c(^|i)>^(fc|0. 

x(k\t),ak\t),5(k\t),z(k\t)J 


k=t-T+ 1 
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Three  Tank  System:  Fault  Detection  -  III 


measured  output 


estimated  output 


•  <j> i  leak  in  tank  1  for  20s  <t<  60s 

•  fa  valve  V\  blocked  for  t  >  40s 

•  -f-  logic  constraint  [hi  <hv]=><f> 2  =  0 
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Research  Program 


•  Modeling  Language  and  Compiler 

—  Express  models  and  specs  with  “standard” 
vocabulary 

—  Translate  into  MILP/MIQP 

•  Verification/Controllability/Reachability 

—  Use  tools  from  Polyhedral  Computation 
—  New  case  studies 

•  State  Estimation/Fault  detection 

—  Dual  of  Model  Predictive  Control 
—  Convergence  properties  of  estimator 

•  Optimization  Algorithms 

—  Sparse  QP  solver 
—  New  Branch-and-Bound  strategies 


•  Model  Reduction/Approximation 
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HYSDEL 

(HYbrid  Systems  DEscription  Language) 


•  Describe  Hybrid  Systems  in  a  compact  way 

—  Propositional  Logic 
—  Dynamics 
—  Constraints 

•  Automatically  generate  MLD  models 

•  MLD  model  is  not  unique  in  terms  of  the 
number  of  auxiliary  variables  optimize 
model:  #  binary  variables  =  min.! 

•  Reduction  of  model  complexity  for  models 
containing  purely  logic  relations 

(Truth  Tables  =>-  Convex  Hull) 


ZBrMi 
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HYSDEL  Example 
Automotive  Active  Leveler 


X  Description  of  variables  and  constants 

state  £,h,xll,xl2; 
input  d,dc,dev; 

const  OTh,  0T1,  ITh,  in; 
const  Ml,H2,H3,M4,mi,i»2,ra3,m4; 
const  e; 

const  Ts ,cbar, evbar , eats ,  cmax , cmia.evmax, evmin; 

X  Variable  types 

real  f  ,li,zt,z2,z3,d,dc,dev; 

logic  dl,d2.d3.d4.d5,d6,d7,d8.d9,dl0.dll.di2.dl3,dl4; 


X  Relations 


dl  =  -Cf-ITh.  <=  0.  Hi,  ml,  e>; 

d2  =  {f-m  <=  0,  M2,  m2,  e>; 

d3  =  -Cf-OTb  >=  0,  M3,  m3,  e>; 

d4  =  -Cf-OTl  >=  0,  H4,  m4,  e>; 

d5  -  "xll  &  d3;  X  Should  be  accepted  also:  d5=(l-xli)ftd3,  d5“(l-xll)*d3 

d6  =  ill  &  "dl; 

d7  -  xl2  ft  d2; 

d8  -  "xl2  ft  "d4; 

d9  =  dS  I  d6; 

dlO  -  d7  I  d8; 

dll  -  "xll  ft  x!2; 

dl2  =  xll  ft  "xl2; 

d!3  =  "("d9  ft  "dlO  ft  di4) ; 

dl4  =  xll  I  xl2; 

zl  -  dll  *  (cbar  +  dc)  {cmax.cmin.e}-; 
z2  =  dl2  *  (evbar  +  dev)  {evniax, evmin,  e}-; 

z3  =  (eats  *  f  +  (1  -  eats)  *  h)*dl3  {10»0Th,10*0Tl,e}; 

X  Other  constraints 

must  xll  +  xl2  <=  1; 
nust  " (d9  ft  dlO) ; 
must  "(dll  ft  dl2); 

X  Update 
update  f  =  z3; 

update  h  =  h  +  Ts  *  (d  +  zl  +  z2) ; 
update  xll  =  d9; 
update  x!2  =  dlO; 
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General  Transformation  of  Truth  Tables 
into  Linear  Integer  Inequalities 


4 

A8<B,  8  €  {0,  l}n 

Theorem: 

The  polytope  P  =  {8  :  A8  <  B}  is  the  convex 
hull  of  the  rows  of  the  truth  table  T 

p  =  conv(T) 

(Bemporad  and  Mignone)  □ 

Every  logic  proposition  can  be  translated  into 
linear  integer  inequalities 
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Truth  Tables  =>  Linear  Integer  Inequalities 


<5i 

S2 

S3 
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0 

Example:  “AND” 
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/  t 


(0,0,0) 
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(1,0,0) 


(1,1,1) 


0(14,0) 


-  8, 


conv ([8],  [i],  [o],  [i])  = 


—Si  +  S3  <  0 
—62  +  S3  <  0 
+  S2  —  ^3  <  1 


Algorithm  to  compute  convex  hull:  lrs 

(David  Avis,  McGill  Univ.)  ftp://mutt.cs.mcgill.cs/pub/c 


Other  algorithms:  qhull,  chD,  Hull,  Porto,  cdd. 
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Research  Program 


•  Modeling  Language  and  Compiler 

—  Express  models  and  specs  with  "standard” 
vocabulary 

—  Translate  into  MILP/MIQP 

•  Verification/Controllability/Reachability 

—  Use  tools  from  Polyhedral  Computation 
—  New  case  studies 

•  State  Estimation/Fault  detection 

—  Dual  of  Model  Predictive  Control 
—  Convergence  properties  of  estimator 

•  Optimization  Algorithms 

—  Sparse  QP  solver 
—  New  Branch-and-Bound  strategies 

•  Model  Reduction/Approximation 
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Formal  Verification  of  Hybrid  Systems 


•  MLD  Model: 

x(t  +  1)  —  Ax(t )  +  B\w(t)  +  B2S(t )  +  B3z(t) 
y(t )  =  Cx(t)  +  Diw(t)  +  D2S(t)  +  D3z(t ) 
E2S(t)  +  E3z(t)  <  E!w(t)  +  E4x(t)  +  E5 


Verification  Problem  VP  # 1 :  Find  max  range  for  y(t) 
Vt  >  0  ,  w(t)  G  W  ,  x(0)  G  Af(0) 

Verification  Problem  VP  #2:  Vw  g  W  and  Vx(0)  G 
Af(0)  verify  that  x(t)  g  Afs  (set  of  safe  states) 


•  Simple  solution  (VP  #1):  Solve  VT  >  0 

max  y(T)  =  Cx(T)  +  jDitu(T)  +  D28(T)  +  L>3*(T) 


subj.  to 


x(0)  G  Af(0) 

w(t)  G  W,  0<t<T 

x(t  +  1)  =  Ax(t)  +  £iw(£)  +  B28{t)  +  B3z(t) 
E25(t )  +  E3z(t)  <  Eiu(t )  +  J&4x(t)  +  £J5 


•  Similarly  for  VP  #2:  feasibility  test  with  x(T)  0  Ai 

•  IMPRACTICAL  !  (even  for  polyhedral  sets  W,  ATSf 
A'(O),  because  of  integer  constraints) 
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Verification  Algorithm 

•  RMK:  when  S(t)  =  const  ,  system  behaves  linearly. 

•  IDEA:  hypothetically  partition  MUc  =  \jf=i  Ci  (N  <  2rt ) 
where  xc(t)  eCi  o  5(t)  =  Si 


subj.  to  < 


v 


x(k  +  1)  =  Ax(k)  +  B\w(Jt)  +  B25(k)  +  Bsz(k) 
E2S(k )  +  £32 (fc)  <  Eiw(k )  +  Et\x(k)  +  E5 

w(k )  ew,  k  =  0,...  ,t 

6(k )  =  Si,  k  =  0, . . .  ,t  -  1,  <5(t)  e  {0,  l}r' 

*(0)  6  *(0) 


•  Reachable  set  implicitly  defined  by  linear  inequalities 
(Hp:  W,  X(0),  Xs  polytopes) 
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Example:  Verification  of  an 
Automotive  Active  Leveler 


System 


'  "f 


compressor 


Escape 
'  ■  Valve 
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•  Controller 


f<oiy ^|| 
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->f=o 

in 
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i 
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f>OTl^ 

y^OTh 
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•  For  a  given  set  of  initial  conditions,  variation  of 
parameters,  and  road  conditions: 

(VP#1)  Find  worst  performance 

|(VP#2)  Verify  that  height  is  within  desired  limits 
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Verification  of  an  Automotive 
Active  Leveler  -  II 


•  System  can  be  transformed  into  MLD  form 


Continuous  states 

2 

Logic  states 

2 

Disturbance  inputs 

3 

Boolean  auxiliary  variables 

14 

Continuous  auxiliary  variables 

3 

•  Results  and  computational  times: 

-44.54  <  h(t)  <  25.00  4  m  I  Pentium  II  400  MHz 

25  m  SPARCstation  20 
(M-code) 


—  No  use  of  fast  iterative  reach-set  projection 
—  Use  of  rectangular  approximation  of  initial  sets 

•  Using  HyTech  (Stauner  eta/.,  1997): 


-47  <  h(t)  <27  62  m  SPARCstation  20 


Analytical  solution  (Elia  and  Brandin,  1998): 


-43  <  h(t)  <  23  | - by-hand 

Note:  exact  limits  in  discrete-time  7^  continuous-time 

mmmmmm  m  1 

-  Mm  W  jgajg  TmohtU* abm  Moa&MAufcj  - — — - 


!  JEMf/tm&fghehm 
T+cluU*ttbm  Hoa h 
ZOrMi 


Research  Program 


Modeling  Language  and  Compiler 

—  Express  models  and  specs  with  “standard” 
vocabulary 

—  Translate  into  MILP/MIQP 

Verification/Controllability/Reachability 

—  Use  tools  from  Polyhedral  Computation 

—  New  case  studies 

State  Estimation/Fault  detection 

—  Dual  of  Model  Predictive  Control 

—  Convergence  properties  of  estimator 

Optimization  Algorithms 

—  Sparse  QP  solver 

—  New  Branch-and-Bound  strategies 

Model  Reduction/Approximation 
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Branch  &i  Bound  for  Optimal  Control 
(and  Fault  Detection)  of  Hybrid  Systems 

(A.Bemporad,  D.  Mignone,  M.Morari) 

Key  Idea: 

1.  Optimal  sequences  {  5(0),  5(1),  ...  ,  5(T)  } 
have  FEW  switches 

*  8(t) 


example:  [5(t)  =  1]  ■*-»  [ x(t )  >  0] 


2.  Order  the  sequences  by  number  of  switches: 

0000,  1111;  0001,  0011;  0010,  0110;  ... 
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Branch  Si  Bound  for  Hybrid  Systems 


•  First  explore  combinations  of  6  with 
low  number  of  switches  (“outside  first” 
strategy) 

•  Assumption:  global  optimum  is  reached 
early  =>  large  subtrees  are  fathomed. 

•  If  not  enough  available  computation  time 
look  only  for  sequences  having  a  number  of 

switches  <  Kmax 


mmmgmm  m  i 
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depth  first 
breadth  first 
outside  first 


Three  Tank  System:  Control  -II 
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Trajectories  using  suboptimal  MIQP  solutions 
Limiting  the  number  of  QPs  at  each  time  step 


sotfnwtt  40 


20  Dm.  30 


global  optimum 


limiting  #  QPs 
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Premise 


•  The  discussed  problems  are  important  but 
inherently  difficult 

•  All  useful  techniques  must  involve  significant 
off-line  and/or  on-line  computation 


Results 

•  New  System  Type: 

Mixed  Logical  Dynamical  (MLD)  System 

•  Many  practical  problems  can  be  represented 
in  MLD  form 

•  Control,  estimation,  and  verification 
require  solution  of  Mixed-Integer  Linear  (or 
Quadratic)  Programs  (MILP,  MIQP)  for 
which  efficient  techniques  are  becoming 
available 
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ROLLOUT  ALGORITHMS: 
PERFORMANCE  ANALYSIS 


Rollout  Algorithms:  Performance  Analysis 
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DYNAMIC  PROGRAMMING  / 
DECISION  AND  CONTROL 
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DP  and  NDP 
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ROLLOUT  POLICIES 
(Tesauro  1996) 


in 


> 

O 
■  ■■■ 

o 

On 

m  CL 

^  <i) 

<5  co 
1  £ 
O  3, 


3** 

"***  "T"!  C 

®  (Q  § 
qO'^ 
O  &  « 
0)  ts  3 

>  O  £ 

<  S  (/) 


T5  “ 

<5  ° 

% 

O  © 

v.  V) 

<D  5 

4-J»  .Q 

m  — 
£  « 
l.i 

2  Q. 
Q.  O 


CO  3  c 
</,</>  O 
CC  (D  "J3 

O  E  £ 

O  CD 

co  ±: 


(o 


Rollout  Algorithms:  Performance  Analysis 


X 

0 

c 

o 

XX 

> 

"55 

o 

CD 

H 

o 

o 

Q. 

Cl 

0 

15 

CD 

CIS 

jQ 

<5 

TS 

CCS 

0 

CD 

C 

CCS 

o 

o 

o 

Id 

V\ 

0 

“SI 

0 

IS 

0 

+-* 

0 

■£H 

£ 

0 

o 

0 

0 

w 

0 

CD 

< 

T3 

-C 

o 

0 

x- 

0 

0 

0 

\ _ 

E 

o 

o 

it 

0 
■  MM 

Li 

0 

3 

o 

0 

c 

SI 

0 

4-» 

s*-^' 

0 

o 

0 

■  ■■ 

.D 

0 
■  ■MB 

JZ 

HH 

3 

0 

■  mm 

-C 

<u 

0 

■M 

0 

JI 

4-» 

0 

c 

3 

>< 

0 

c 

"D 

0 

c 

0 

sz 

■i"- 

■> 

0 

0 

o 

0 

0 

0 

0 

(/) 

#  •  • 


Rollout  Algorithms:  Performance  Analysis 


> 

a  co 

Op 

tTLU 

3  Q_ 

Qo 

cc 

O0- 

oc 


H  to 

I  ■o 

g  0 

O  1L 

a 
co 

(D 
a*. 

.2 
O 

0 

3 
0 


0 
£ 

3 

O 

.2 

3 
0 
£ 

0  -C 

£  ° 

O)  ^ 
-£ 

-x  .£ 

O  -H 
o  o 
—  0 

TJ  £ 

}=  o 
§  o 

C  2= 
O  0 
UL  C/) 


O) 

£ 

IM 

£ 

£ 

0 

a 

0 


0 

0 

3 

0 


a  £ 
0  o 

%  o 

£  o 

"4=  XJ 
0 

O  « 

0  5 
0  0) 


>*  0 
.2  -* 

o  (/) 
a-£ 

ti  0 

oSP 
=  0  o> 
0.2^ 


22  CC 


.12  - 

jj  £  0 

agS 

8^.J2 

c  <D  i= 

+-  »  s 

c  to  ® 
a>  J3  X 

E  On 

or  „ 

>  +-  3 

o  c  > 

o.°K 

£  8  » 
>  2 


CO  o 


>QZ 
w  ... 


o  £  0 
=  as 
o  c  0 

CL  JE  I— 


U> 

0 

■N 

0 

a 

£ 

0 
■  ■■ 

0 

a. 

x 

0  a) 

>  E 
.£  0 

-*-»  — 

0  -9 
o  2 
a.  q. 

75  g) 

£  = 
0 

0  d) 
X2  jr 

3  a 
C/)  0 


(0 

32 

w 

•c-« 

'</> 

I- 

=  1 

O  g 

«  » 
(/)  ® 
(S  ^ 
o  a 
w  <0 
3  t 
c  o 

(0  CO 
CO  - 

2  3 

CD  jc 
(0  O 
■E  c 

CD  O 
CD  D 


Rollout  Algorithms:  Performance  Analysis 


PERFORMANCE  QUESTIONS 


c^- 

o 

0 
■  warn 

3 

0 

0 

(0 

cd 

n 

0 


CD 

0 

£ 

E 

3m 

■  mmm 

CD 

Q. 

JJ 

O 

0 

-Q 

O 

-C 

0 

O 

CO 

3 

O 

E 

O 

5 

5 

o 

o 

i 

I 

• 

• 

0» 

u 

cd 

0 

03 

O 

o 

o 

CD 

N 

'</> 


O 

ol 


r* 

n^rffr^ 

MQ 

C 

o 


CO 


CD 

0) 

CD 

JC 


(0 
K  SIMM  MWW 

m  ® 

O  (0 

■  man  m m 
■I mdt  mam 

S2 ‘55 

$**>  HM 

fefW  JMX 

CD  5 

f  o 

******  jl  Ik 

«$*«*  w 

CD  «*. 

Emm 

O 

U>JC 

c  g 

«WMtt 

O  O 
13  (0 
h—  CD 

o  o  ,ti 


rr^m^ 

SWM 

CD 

jQ 

H  UMAX 


«8  £ 

a  cd 

CD  CD  C 


L- 

CD 


t.'SSSSRJ'.'SwSSJ'JSKSSS 


Rollout  Algorithms:  Performance  Analysis 


THE  BREAKTHROUGH 
PROBLEM 
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AVERAGE  PERFORMANCE 

ANALYSIS 


•  •  • 
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OPTIMAL  BREAKTHROUGH 

PROBABILITY 
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GREEDY  BREAKTHROUGH 

PROBABILITY 
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ROLLOUT  BREAKTHROUGH 

PROBABILITY 
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ROLLOUT  TO  GREEDY  RATIO 
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ROLLOUT  WITH  LIMITED 
LOOKAHEAD  m 
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EXPERIMENTAL  OBSERVATION 
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ROLLOUT  TO  GREEDY  RATIO 
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The  Behavioral  Approach 

to 

Systems  Modelling 
and 

Control 

Jan  C.  Willems 
University  of  Groningen 
The  Netherlands 
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FEEDBACK  PASSIVATION  DESIGNS 
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ACTIVATION  OF  CONCEPTS 


Over  the  last  ten  years  several  descriptive 
concepts  have  been  activated  into  constructive 
design  tools  such  as 


•  Lyapunov  functions  — >  Control  Lyapunov 

functions  (CLF) 

•  Total  stability  Input/state  stability  (ISS) 

•  Passivity  -4  Passivation 

•  Optimality  — >  Inverse  optimality 

leading  to  design  procedures  like 

•  Cascade  passivation 

•  Backstepping:  robust,  adaptive,  stochastic 

•  Forwarding 

•  Passivation  using  LMFs 


FOUNDATIONS: 

absolute  stability:  Lurie  (1951), 

passivity:  Popov  (1960), 

PR  Lemma:  Yakubovich  (1962),  Kalman  (1963), 

passivity  and  small  gain:  Zames  (1966), 

inverse  optimality:  Kalman  (1964),  Moylan 
and  Anderson  (1973), 

dissipative  systems:  Willems  (1972),  Hill  and 
Moylan  (1976), 

conicity  and  margins:  Safonov  and  Athans 
(1977),  Safonov  (1980). 

DESIGN  TOOLS  AND  PROCEDURES: 

Lyapunov  design:  Meilakhs  (1976,1978), 
Jurdjevic  and  Quinn  (1978),  Corless  and 
Leitmann  (1981),  Qu  (1998), 

control  Lyapunov  functions:  Artstein  (1983), 
Sontag  (1989), 


backstepping:  Tsinias  (1989),  Byrnes  and 
Isidori  (1989), 

ISS  and  ISS  small  gain:  Sontag  (1989),  Jiang, 
Teel  and  Praly  (1994),  Teel  (1996), 

passivation:  PK  and  Sussmann  (1989),  Byrnes, 
Isidori  and  Willems  (1991), 

adaptive  backstepping:  Kanellakopoulos,  PK 
and  Morse  (1991),  Krstic,  Kanellakopoulos  and 
PK  (1992), 

robust  backstepping:  Freeman  and  PK  (1992), 
Marino  and  Tomei  (1993), 

ISS  and  Hoo  backstepping:  Jiang,  Teel  and 
Praly  (1994),  Pan  and  Ba§ar  (1998),  Krstic  and 
Li  (1998), 

stochastic  backstepping:  Krstic  and  Deng 
(1998),  Pan  and  Ba§ar  (1999), 

forwarding:  Teel  (1992),  Mazenc  and  Praly 
(1996),  Jankovic,  Sepulchre  and  PK  (1996). 


In  this  lecture 


•  Locally  optimal  backstepping  design 

Ezal,  Pan  and  PK  (1997) 

•  Passivation  redesign 

Arcak,  Seron,  Braslavsky  and  PK  (1998) 

•  Observer-based  feedback  passivation 

Larsen  and  PK  (1999) 

•  Passive  observer  design 
Arcak  and  PK  (1999) 

•  Example:  Ship  steering 

Arcak,  Fossen  and  PK  (1999) 


CLF  Construction 


Backstepping  Lemma  If  V (X)  is  a  CLF  for 

X  =  F(X)  +  G(X)u 

and 

L(F+Ga)V  =  LfV(X)  +  LGV(X)  a(X)  <  -ff(X) 
then, 

V+(X,x)  =  V(X)  +  (x-a(X))2 
is  a  CLF  for 


X  =  F(X)  +  G(X)x 
x  =  f(X,x)  +  g(X,x)u  g(X,x)^0. 


Example 


X  =  X3  +  x 

X  =  U. 


V(X)  =  X2,  a(X)  =  —X3  —  X 
V+(X,  x)  =  X2  +  (x  +  X  +  X3)2. 


Xo  —  l'o(Xo)  +  Go(Xo)xi 
XT  =  fi  (X,xt  )  +  gi  (X,xt  )x2 

X2  =  f2(X)X1,X2)  +  g2(X)Xl,X2)X3  (SF) 

•  • 

•  _____  • 

xn  fnfX^x)  ~\~  Qn (X, x)tl . 

9i(X,  •  •  •  ,Xi)  ^0 
L(F0  +  G0oc0)^o(Xo)  <  — Co(Xo) 

Step  1.  Set  X  =  X0,  x  =  xi ,  V(X)  =  V0(X0),  apply 
backstepping  lemma:  V+  =:  Vi  (X,  xt  ) . 

Step  i.  (i  =  2,  •  •  •  ,  n) 


X  = 


Xo 


Xi-T 


x  —  Xi  — >  V+  — :  Vt(X,xt  ,  •  •  •  ,xt). 


Step  n  — > 


Vrt(X,  X,  ,  •  •  •  ,Xn) 


Locally  Optimal  Backstepping  Design 

K.  Ezal  et  al.  -  CDC  97  Best  Student  Paper 


x  =  Ax  +  }{x)  +  G  i{x)w  +  G2(x)u  (strict-feedback  class) 

Robust  backstepping: 


Desire  Local  Tioo  -Optimality:  x  =  Ax  -f  B\w  +  B2U 


Requires  r(0)  =  R ,  qXx{ 0)  =  Q,  Bi  —  Gi(0) 

^>0,  Q  =  >0-a  priori  specified 


Solution:  Cholesky  factorization 

PA  +  A'P  +  P  {'r~3B1B,1  -  B2R-lB'2 )  P  +  Q  =  0 


-61 

O' 

P  —  V 

_o 

<5n. 

L, 

L  = 

1  0 

■aii  1 


:  '  •  •  ' '  •  0 

-Gfn  — 1,1  ~~an—  l,n—  1  b 


Result:  Locally  Uoo  -optimal,  globally  inverse  optimal 


Stability  region  of  an  Hoo  design 

*&1  =  X2  -f  X2  "}~  U) 

X2  =  U 

Linearized  system: 

%1  =  X2  -f  W 

X2  =  U 

Cost  functional: 

rOO 

J  =  J  [®?  +  x\  +  u2  —  72tn2]  dt 
Desired  attenuation:  7  =  5  >  7*  =  1.27 


Optimal  linear  control: 

u  =  -B'Px  =  -1.06xi  -  1.78x2 


Passivation  Redesign 

Arcak,  Seron,  Braslavsky  and  PK  (1998) 


Ho  Nominal  plant 


x  =  ®(x)  +  r(x)v 


Input  unmodeled  dynamics 


t  =  fA(f,,u) 
v  =  h.A(i,u). 


A  is  0-GAS  and  3  pos.def.  TT ( £,)  : 
ri(y  ^  — Au2  T  vu,  A  >  0  (excess  of  passivity) 


This  includes  stable  linear  A (s)  with 


Re{A(jeu)}  >  A. 


Find  pos.def.  V(x)  such  that 

UV(x)  <  (l-rV(x))2  Vx^  O. 


Then, 


=»  GAS, 


Proof:  Shortage  vs.  excess  of  passivity 


Shortage  of  passivity  1  /k: 
kV  =  kLcpV  +  kLrVv  <  k(LrV)2 +  kLrVv 

"uv 


Excess  of  passivity  A:  fl  <  — Au2  +  vu 

U(x,  0  =  TcV(x)  +  TT(£,)  =>  U  <  0  Vy  ^  0. 

X  =  0=^u  =  0=^£,  =  0is  the  largest  invariant 
set  by  0-GAS  of  A. 


When  does  such  V(x)  exist? 
Iff  3  CLF  V(x)  such  that 


lim  sup 
x-*o 


UV(X) 

(LrV(x))2 


<  +oo. 


Construction  of  V(x): 

Find  positive  0(-)  such  that 

(i)  uvfxi  •  <  9(V(x)) , 

W  (LrV(X))2 

(ii)  limt-ioo  Jo  6(s)ds  =  +°°  . 
Then,  |v(x)  =  fjf(x)0(s)ds 


(L) 


Lemma  When  Ho  is  strict  feedback  (SF),  our 
backstepping  CLF  Vn(X,xi ,  ••  •  ,xn)  constructed 
from  Vo(Xq)  with 


Lf  o  +  G  0  OCq  Vo(X0)  <  -ffo(Xo) 
satisfies  (L)  if  Vq  and  Oq  have  pos.def.  Hessians. 


Passivation  Redesign  Example 


X  =  X3  +  X 


X 

V 


V 


12(s  +  35)(s  +  20) 
7(s  +  40)(s  +  30) 


u  =  A(s)u. 


Excess  of  passivity  Re{A(jcu)}  >  A  =  1 . 


V0(X)  =  X2,  oc0(X)  =  —X  —  X3,  a0(X)=2X2 


3 \2 


Vi(X,x)  =X2  +  (x  +  X  +  X5) 


L©Vi(X,x)  3  9  2 


(LrVi  (X,x)) 


a<i  +  4Vl(X>x) 


Choose  9(Vi )  =  §  +  fVf ,  k— ^  — 1, 


u  =  -kLrV  =  -0(Vi  )LrVi  =  -20(Vi  )(x+X+X3). 


Feedback  Passivation  (FPR)  Design 

In  design  we  use  state  feedback 

u(x)  —  oc[x)  +  |3(x)u 


to  make  H  passive  for  the  input-output  pair 
(u,y): 


This  is  possible  if 

•  H  is  relative  degree  one: 

u  appears  in  the  q-equation, 

•  H  is  (weakly)  minimum  phase: 

the  “zero  dynamics”  which  remain  when  the 
output  is  kept  at  zero,  q(t)  =  0,  are  stable. 


Observer-based  FPR  Design 

Larsen  and  PK  (1999) 

x  =  fo(x,0)  +  g(x,-y2)yi 
£,  =  A £,  -P  Bu 
yi  =  C-j£,  V2  =  C2£, 

•  Design  an  asymptotically  stable  observer 

t  =  At  +  Bu  +  L(y2  -#2) 

$2  =  C2t 

AjM  +  MAo  =  -I,  A0  :=  A  +  LC2,  M  >  0. 

•  Design  a  state  feedback  gain  K  to  stabilize  the 
linear  part  of  the  system  with  Ak  :=  A  —  BK 
Hurwitz  and  (A^,  B,  C-| )  passive  by  solving 

A^P  +  PAk  =  -Q  <  0  P  >  0 
PB  =  C}. 

•  Design  a  stabilizing  u  from  the  CLF 

W(x,  e)  =  V(x)  +  +  — eTMe. 
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follows  from  the  LaSalle  Invariance  Principle  using  the  detectability  of  (Ak,  Q )• 


State  feedback  FPR  design 


Observer  feedback  FPR  design 


Passive  Observer  Design 

Arcak  and  PK  (1999) 

x  =  Ax-ip(x)  +  T(y>u) ,  y  =  Cx, 

Observer: 

it  =  a£--4>(£)  +  "yCv,^)  +  L(y  -  at) 


Observer  gain  matrix  L  to  be  designed  for 


iMx) 


ipiUi) 
4*2  (X2) 


n  (Xn) 


(ND)  =*>  a  [^i(xi)  -  rki(xi  —  cr)]  >  0  ( sector ) 


Error  dynamics: 

e  =  (A  —  LC)e  —  (ip(x)  —  4>(£)) 

Indicator  matrix  K  locates  the  nonlinearities: 

K  =  diag(ki,  •  •  •  ,kn) 

f  0  if  iIh(xi)  =  0 
kt  =  < 

1  otherwise 

Note:  PK  =  K  does  not  imply  P  =  I. 

z  :=  Ke 

<t>(t,z)  :=rMx)  -xKk) 

zT4)(t,z)  >  0  (sector) 

e  =  (A  —  LC)e  —  Kcf>  (t,  z) 
z  =  Ke . 


Observer  Design 


V 

3- 

A 

e  =  (A-LC)e  +  Kv 

z  =  Ke 

Z 

4>(t, -) 

Find  P  =  PT  >  0  and  L  such  that 
(A-LC)tP+P(A-LC)  <  0  and  PK  =  K  . 

Then,  whenever  e(t)  exists,  it  satisfies 

|e(t)|  <  k|e(0)| . 


Proof:  PR  and  sector  properties. 


Control  Law  Design 

x  =  Ax  — iMx)  +y(y,u) ,  y  =  Cx, 
it  =  +y(-y,u)  +  L(y  -  Cit) . 

Find  u  =  oc{y ,  it)  to  guarantee 

|x(t)|  <  max  {  (3(|x(0)|,t)>  9  (supT6[0>t)|e(T)|)  } .  (ISS) 

Main  Result 

If  (PR)  and  (ISS)  a, re  satisfied,  then  (CL)  is  GS. 

If,  in  addition,  (K,A  — LC)  is  detectable  and, 

\(Vs  &re  either  strictly  increasing  or  identically 
zero,  then  (CL)  is  GAS. 

Proof: 

(PR)  implies  |e(t)|  <  k|e(0)|  and,  with  (ISS), 
proves  GS.  From  (PR)  and  Barbalat’s  lemma 

z(t)  =  Ke(t)  — )  0. 

(CL)  is  autonomous  and  z  =  0  dynamics  are 
e  =  (A  -  LC)e.  Detectability  implies  e(t)  0 
and  GAS  follows  from  ISS. 


Feasibility 

(PR)  is  an  LMI  in  LTP  and  P  =  PT  >  0: 

(A  -  LC)tP  +  P(A  -  LC)  PK  -  K 
KP-K  0 

Feasibility  depends  on  K,  the  location  of  ipi(xt). 

Example  y  =  xi 

Xt  =  X2 

x2  =  — 4>2(x2)  +y[v)  +  u* 

0  0 

ipi  (xt  )  =  0  — >  K  = 

0  1 

(PR)  satisfied  with  P  =  I,  L  =  [  1  1  ]. 

Evaluate  V  for  V  =  \  ef  +  \  e2  along 
ei  —  — ei  +  c2 

e2  =  -d  -  (^2(^2) -^2(^2)). 

V  =  -e,  -  e2(rjj2(x2)  -4>2(x2  -  e2))  <  -ef. 


Example:  Ship  Steering 

Arcak,  Fossen  and  PK  (1999) 

ii  =  J(r|)v 

Mv  =  —  D('v),v  —  g  (rj)  +  t. 

Si('Vi) 

D(v)v  =  Dov  +  j 

M^) 

Do  =  Dj  >  0  ,  5i(-)’s  are  ND. 

Observer: 

f|  =  Kn^  +  h-ft), 

—  —  D(^)^  —  g(rj)  +  t+  JT(ri)(r|  —  fj) 


V(il,v)  =  ]jf\Tf\  +  Vmy  , 
-^PTfi  -  ^•vtD0v-vt[6(v)  -  6(v 


<  -ah  d 


vDoi . 


V 


Observer  Based  Control  Design 


V  =  *1 

X-\  —  X2  ~h  Xj 

X-2  =  X3  —  X2(l  +x|) +U 

X3  =  2u . 


'  0 

1 

0  " 

A  = 

0 

-1 

1 

C  =  [1 

1  0  0] 

0 

0 

0 

0 

y2 

— 

*2 

y(y,u)  = 

u 

0 

2u 

->  K  = 


0  0  0 
0  1  0 
0  0  0 


(PR)  satisfied  with  L  =  [  3  3  1  ]T. 


ISS  Control  Law  Design  by  Backstepping: 


V  =  *1 

Xi  —  X2  +  xf 

x2  =  X3  —  X2(1  +  x2)  +U 

x3  —  2u . 

Relative  degree  two  +  ISS  zero  dynamics 
— >  two  steps  of  observer  backstepping: 

Step  1.  xi  =  x2  + 

cxi  (x-| )  =  — xf  —  xi  —  x-\  ,  ui  >  0. 

XT  =  — Xi  —  TLl  Xl  +  e2  +  (&2  —  <*1  ) 


Step  2.  Z2  :=  £2  —  (*1 )  =  X2  —  oc-i  (xi )  +  ez 

,  /\  \  3ai 

z2  =u  +  q(xi,x)  -  - — e2 

0x1 


Proof  of  ISS: 


ISS  Lyapunov  function:  V(xi ,  zz)  +  z\ 
— »  (xi ,  X2 )  subsystem  is  ISS  with  input  e2- 


ISS  of  zero  dynamics: 

r|  :=  x3  -  2x2 

f|  =  — 2rj  —  2x2  +  2x2  (2D) 

Cascade  of  ISS  systems  is  ISS: 


ISS  ISS 
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•  Polynomial  splines,  in  particular  cubic  splines,  are  the  most  well  known. 
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Configuration  spaces  for  mechanical  systems  are  usually  non-Euclidean 


•  Notations: 

For  a  real-valued  functions  s,  defined  on  the  time  interval  [0,T], 


6 


(«fe,  Pi,oi  Pi,ni  Vi,  k,  are  given  real  constants) 


Theorem  [Ahlberg  Nilson  and  Walsh,  1967] 

Given  the  differential  operator  L  and  the  partition  A  of  the  time  interval  [0,  T], 
there  exists  a  unique  generalized  spline,  for  each  set  of  boundary  and  interpolation 


To  find  this  unique  spline  it  is  enough  to  determine  2 mn  unknowns,  correspond- 


An  optimal  control  problem  with  interpolation  conditions 

' 

distinct  points  x0,xt  in  the  state  space; 
n  —  1  real  constants  a*.,  k  =  1,  2, . . . ,  n  —  1; 


(Martin,  Enqvist,  Tomlinson,  Zhang  -  1995),  ( Rodrigues ,  Silva  Leite,  Simoes  -  1999 ) 


•  The  output  function  corresponding  to  the  previous  optimal  control  problem  is  a 
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(A)  =  {0}  cr(A)  =  {1,-10} 


•The  optimal  control  problem  with  arbitrary  output  function  is  similarly  related 
to  another  class  of  splines,  the  L-splines.  This  case  allows  uneven  smoothness 


•  Cubic  spline  in  JRm  -  a  variational  approach 
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The  Euler-Lagrange  equation  associated  with  this  problem  is: 


(double-bracket  equations) 


Difficulties: 

•  It  remains  an  open  problem  to  describe  solutions  of  these  double-brackett  equa- 
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•  Where  are  the  connections  with  optimal  control? 

t  It  is  possible  to  formulate  a  number  of  optimal  control  problems  which  are  equiv¬ 
alent  to  the  variational  problems  associated  with  non-Euclidean  splines. 
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Case  III  -  nonholonomic  constraint,  with  drift 


•  The  De  Casteljau  algorithm  -  a  geometric  construction  of  splines 

(De  Casteljau,  Technical  Report,  Citroen/Paris  -1959) 

•  The  classical  De  Casteljau  algorithm  is  a  geometric  construction,  whereby  two 


Illustration  for  tlie  construction  of  tli©  cubic  polynomial. 
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Relationship  between  control  points  and  initial/final  velocity/acceleration: 
p3(0)  =  Xq ,  P3(l)  =  #3? 

p3(0)  =  3(a?i  —  aj0)?  P&0-)  —  —  ^2)5 


10 

05 


•  The  power  of  this  algorithm  lies  in  the  fact  that  it  can  be  easily  generalized  from 
JRm  to  other  spaces,  as  long  as  the  linear  interpolation  process  is  suitably  redefined. 


CO 

C4 


•  Generalized  De  Casteljau  algorithm  on  Lie  groups 


Relationship  between  control  points  and  initial/final  velocity/cov.  acceleration: 
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•  Main  difficulties  in  implementing  the  De  Casteljau  algorithm: 
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•  So,  for  SO (3),  the  De  Casteljau  algorithm  to  produce  polynomial  curves  and 
spline  curves  is  easy  to  implement. 
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Generalized  De  Casteljau  algorithm  on  spheres  S 


So,  the  De  Casteljau  algorithm  on  spheres  can  be  easily  implemented. 
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•  The  class  of  interpolating  curves  defined,  by  a  variational  principle  or  by  a  Hamil¬ 
tonian  approach  does  not  seem  so  computationally  tractable  as  those  developed  by 
the  De  Casteljau  algorithm. 
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Good  reasons  to  study  the  geometry  of  splines: 
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